/*
 * yytSGA.cpp
 *
 *  Created on: 10-04-2011
 *      Author: mith
 */

#include<cstdlib>
#include<cstdio>
#include<math.h>
#include<cstring>
#include<time.h>
#include<gsl/gsl_vector.h> // GSL LIBRARY - remember use flag -lgsl -lgslcblas
#include<gsl/gsl_multiroots.h>
#include<gsl/gsl_errno.h> // to interpolation
#include<gsl/gsl_spline.h>
#include"yytSGA.h" // function definitions here!
using namespace std;

// DEFINITIONS
#define MAX_ITERATION 200
#define PRECISION 1e-5
#define NO_OF_EQUATIONS 9
#define DEBUG 1
#define DEEP_DEBUG 0

/*
 *********************************************************
 ************************* GSL ***************************
 *********************************************************
*/

int set (const gsl_vector * x, void *params, gsl_vector * f) {
	parameters p = *((struct parameters *) params);

	variables v;
	v.BDelta = gsl_vector_get (x, 0);
	v.chi = gsl_vector_get (x, 1);
	v.m = gsl_vector_get (x, 2);
	v.d = gsl_vector_get (x, 3);
	v.mu = gsl_vector_get (x, 4);
	v.l_BDelta = gsl_vector_get (x, 5);
	v.l_chi = gsl_vector_get (x, 6);
	v.l_m = gsl_vector_get (x, 7);
	v.l_n = gsl_vector_get (x, 8);

	if(DEEP_DEBUG) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7), gsl_vector_get (x, 8));
	}

	Equa E(p, v);
	variables w;
	E.run(w);

	gsl_vector_set (f, 0, w.BDelta);
	gsl_vector_set (f, 1, w.chi);
	gsl_vector_set (f, 2, w.m);
	gsl_vector_set (f, 3, w.d);
	gsl_vector_set (f, 4, w.mu);
	gsl_vector_set (f, 5, w.l_BDelta);
	gsl_vector_set (f, 6, w.l_chi);
	gsl_vector_set (f, 7, w.l_m);
	gsl_vector_set (f, 8, w.l_n);

	if(DEEP_DEBUG) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5),
				gsl_vector_get (x, 6), gsl_vector_get (x, 7), gsl_vector_get (x, 8));
		printf("###  %24.15e# %24.15e# %24.15e# %24.15e# %24.15e %24.15e# %24.15e# %24.15e# %24.15e#\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			gsl_vector_get (f, 6), gsl_vector_get (f, 7), gsl_vector_get (f, 8));
		printf("\n");
	}

	return GSL_SUCCESS;
}

void print_label(parameters &p) {
	printf("\nsdelta=%f, J=%f, t=%f, U=%f, L=%u\n", p.sdelta, p.J, p.t, p.U, p.L);
	printf("iter \tBDelta \t\tchi \t\tm \t\td \t\tmu \t\tl_tBDelta \tl_tchi \t\tl_m \t\tl_n\n");
}

void print_state (size_t iter, gsl_multiroot_fsolver * s)
{
	printf("%3u %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d) %10.3e (%2d)\n",
			iter, gsl_vector_get (s->x, 0), (int)(log(abs(gsl_vector_get (s->f, 0)))/log(10)),
			gsl_vector_get (s->x, 1), (int)(log(abs(gsl_vector_get (s->f, 1)))/log(10)),
			gsl_vector_get (s->x, 2), (int)(log(abs(gsl_vector_get (s->f, 2)))/log(10)),
			gsl_vector_get (s->x, 3), (int)(log(abs(gsl_vector_get (s->f, 3)))/log(10)),
			gsl_vector_get (s->x, 4), (int)(log(abs(gsl_vector_get (s->f, 4)))/log(10)),
			gsl_vector_get (s->x, 5), (int)(log(abs(gsl_vector_get (s->f, 5)))/log(10)),
			gsl_vector_get (s->x, 6), (int)(log(abs(gsl_vector_get (s->f, 6)))/log(10)),
			gsl_vector_get (s->x, 7), (int)(log(abs(gsl_vector_get (s->f, 7)))/log(10)),
			gsl_vector_get (s->x, 8), (int)(log(abs(gsl_vector_get (s->f, 8)))/log(10)));
}

void print(FILE * stream, double x, variables &v) {
	fprintf(stream, "%f %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e\n",
			(float)x, v.BDelta, v.chi, v.m, v.d, v.mu, v.l_BDelta, v.l_chi, v.l_m, v.l_n);
	fflush (stream);
}

int aPoint(parameters &p, variables &v) {
	const gsl_multiroot_fsolver_type *T;
	gsl_multiroot_fsolver *s;
	int status;
	size_t iter = 0;
	const size_t n = NO_OF_EQUATIONS;
	gsl_multiroot_function f = {&set, n, &p};

	double x_init[n];
	x_init[0] = v.BDelta;
	x_init[1] = v.chi;
	x_init[2] = v.m;
	x_init[3] = v.d;
	x_init[4] = v.mu;
	x_init[5] = v.l_BDelta;
	x_init[6] = v.l_chi;
	x_init[7] = v.l_m;
	x_init[8] = v.l_n;

	gsl_vector *x = gsl_vector_alloc (n);
	for(size_t i=0; i<n; i++) { gsl_vector_set (x, i, x_init[i]); }

	T = gsl_multiroot_fsolver_hybrids;
	s = gsl_multiroot_fsolver_alloc (T, n);
	gsl_multiroot_fsolver_set (s, &f, x);

	if(DEBUG) {
		print_label(p);
		print_state (iter, s);
	}

	do {
		iter++;
		status = gsl_multiroot_fsolver_iterate (s);
		if(DEBUG) print_state (iter, s);
	    if(status) break;	// check if solver is stuck
	    status = gsl_multiroot_test_residual (s->f, PRECISION);
	}
	while (status == GSL_CONTINUE && iter < MAX_ITERATION);

	v.BDelta = gsl_vector_get (s->x, 0);
	v.chi = gsl_vector_get (s->x, 1);
	v.m = gsl_vector_get (s->x, 2);
	v.d = gsl_vector_get (s->x, 3);
	v.mu = gsl_vector_get (s->x, 4);
	v.l_BDelta = gsl_vector_get (s->x, 5);
	v.l_chi = gsl_vector_get (s->x, 6);
	v.l_m = gsl_vector_get (s->x, 7);
	v.l_n = gsl_vector_get (s->x, 8);

	printf ("status = %s\n", gsl_strerror(status));
	gsl_multiroot_fsolver_free (s);
	gsl_vector_free (x);
	return status;
}

/*
 *********************************************************
 ******************** INTERPOLATION **********************
 *********************************************************
*/

//double extrapolation(double &x, double &y, double point) {
//	gsl_interp_accel *acc = gsl_interp_accel_alloc ();
//	gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, POINTS_TO_INTERP);
//
//	gsl_spline_init (spline, x, y, POINTS_TO_INTERP);
//	double w = gsl_spline_eval (spline, point, acc);
//
//	gsl_spline_free (spline);
//	gsl_interp_accel_free (acc);
//	return w;
//}

/*
 *********************************************************
 ************************ MAIN ***************************
 *********************************************************
*/

#define POINTS_TO_INTERP 5
#define RESOLUTION 50.
#define START 20
#define STOP 2
#define OUTPUT_FILE "../mros4_plots/test/test.dat"
#define frand() ((double)rand()/(RAND_MAX+1.0))

int yytSGA (void) {
	time_t t;
	srand((unsigned)time(&t));

	struct parameters p;
	p.J = 1./3.;
	p.t = 1.;
	p.L = 128;
	p.beta = 500;

	FILE * pFile;
	pFile = fopen(OUTPUT_FILE, "w");
    struct variables v, old;
	while(1) {
		p.U = frand()*15.+5.;		// od 5 do 20
		p.sdelta = frand()*0.1; 	// od 0 do 0.2
		v.BDelta = frand()*0.3; 	// od 0 do 0.3
		v.chi = frand()*0.2 + 0.2;	// od 0.2 do 0.4
		v.m = frand()*0.3;			// od 0 do 0.3
		v.d = frand()*0.1;			// od 0 do 0.1
		v.mu = frand()*6. - 3.;		// od -3 do 3
		v.l_BDelta = frand()*10. - 5.;	// od -5 do 5
		v.l_chi = frand()*10. - 5.;	// od -5 do 5
		v.l_m = frand()*10. - 5.;	// od -5 do 5
		v.l_n = frand()*10. - 5.;	// od -5 do 5

		copyTo(v, old);

	    int status = aPoint(p, v); // aPoint take initiate point from v vector and put results back to v
        if(status == GSL_SUCCESS) {
        	if(v.m > 0.05) fprintf(pFile, "EUREKA\n");
        	fprintf(pFile, "# dla: J=%e, t=%e, sdelta=%e, U=%e, L=%d, beta=%f\n",
        	        			p.J, p.t, p.sdelta, p.U, p.L, p.beta);
        	fprintf(pFile, "#\t BDelta \t chi \t\t m \t\t d \t\t mu \t\t l_BDelta \t l_chi \t l_m \t l_n\n");
        	fprintf(pFile, "start: %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
        			old.BDelta, old.chi, old.m, old.d, old.mu, old.l_BDelta, old.l_chi, old.l_m, old.l_n);
        	fprintf(pFile, "wynik: %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
        			v.BDelta, v.chi, v.m, v.d, v.mu, v.l_BDelta, v.l_chi, v.l_m, v.l_n);
        	fprintf(pFile, "\n");
        	fflush (pFile);
        }

	}

	fclose (pFile);
	return 0;
}
